% run this file to recreate Table 5 ("Impacts of Learning and Storage on Macroeconomic Volatility")
% Note that this requires repeated model simulations under different values of parameters/assumptions
% It might therefore take a while to run

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Panel A of the table
% Panel A of the table conditions model simulations on the Kalman filter gain
% The column "misperceiving the Late 2000s", for example, computes simulated moments according to the Late 2000s parameter values
% and then compares these moments with those that would occur if the shocks are from late 2000s but the gain had not changed from the one from the early 2000s, ie. if agents
% had not learned about the change in parameters. As such, the panel seeks to identify the separate impact of learning on model simulations
    
%initialize panel A of the table
table5panelA = zeros(2,2);
comp = [0 1];
track = 0;

% Set up the Kalman filter gain according to early 2000s, for use later

% Parameter Values of oil supply shocks, early 2000s
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p =  0.0655*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8158;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks

% Parameter Values of TFP shocks "early 2000s"
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 0.2707*0.3694*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8466;   % serial correlation of transitory TFP shocks
sigz_t = 0.3694*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


% Call the program that computes the s.s. and the linearization
sstate

% keep filter gain for further use 
gainkeep = P*H/(H'*P*H);

%%%%%%%%%%%%%%%%%%%%%
% Now ready to simulate data

% Set up parameter values related to late 2000s case

% Parameter Values of oil supply shocks, *late 2000s* values
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p = 0.2408*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8294;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;     % std. deviation of innovations to those shocks
% Parameter Values of TFP shocks, *late 2000s* values
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 2.4926*0.2592*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8463;   % serial correlation of transitory TFP shocks
sigz_t = 0.2592*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

% Call the program that computes the s.s. and the linearization
sstate

% keep filter gain for further use 
gaintemp= P*H/(H'*P*H);

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel A of the table:')
disp('Simulating 500 data samples for the *late 2000s parameters - late 2000s gain* solution:')
disp('Will store std deviation(GDP) and std (futures prices) to construct the table :')
simule

%keep the std dev of output and the std dev of futures prices for the table :
table5panelA(1,1) = mvar(iy);
table5panelA(2,1) = mvar(ipof4);
table5panelA


answlearn = 0;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);


disp('Simulating 500 data samples for the *late 2000s parameters - early 2000s gain* solution:')
disp('Will select std deviation(GDP) and std (futures prices) and express them as a ratio of the previous ones to construct the table :')
simule

table5panelA(1,1) = mvar(iy)/table5panelA(1,1);
table5panelA(2,1) = mvar(ipof4)/table5panelA(2,1);
table5panelA

% Now simulate model with "late 2010s" parameters

% Parameter Values of oil supply shocks, "mid 2010s" values
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p = 0.1344*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.6899;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks
% Parameter Values of TFP shocks, "mid 2010s" values
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 0.5624*0.6159*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.9322;   % serial correlation of transitory TFP shocks
sigz_t = 0.6159*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Simulating 500 data samples for the *mid 2010s parameters - mid 2010s gain* solution:')
disp('Will store std deviation(GDP) and std (futures prices) to construct the table :')
simule

%keep the std dev of output and the std dev of futures prices for the table :
table5panelA(1,2) = mvar(iy);
table5panelA(2,2) = mvar(ipof4);
table5panelA


whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 0;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep = gaintemp;
% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Simulating 500 data samples for the *mid 2010s parameters - late 2000s gain* solution:')
disp('Will select std deviation(GDP) and std (futures prices) and express them as a ratio of the previous ones to construct the table :')
simule


table5panelA(1,2) = mvar(iy)/table5panelA(1,2);
table5panelA(2,2) = mvar(ipof4)/table5panelA(2,2);

% OK Done 



%%%%%%%%%%%%%%%%%%%%%%%%%%
% Panel B of the table
%%%%%%%%%%%%%%%%%%%%%%%

% Panel B of the table conditions model simulations on whether inventories are present or not
% The column "early 2000s" for example, computes simulated moments according to the early 2000s parameter values
% and then compares these moments with those that would occur in an economy with the same shock calibration but where inventories are shut off


%initialize panel B Table 5
table5panelB = zeros(2,3);


comp = [0 1];
track = 0;

% First column of the table : *early 2000s* shocks 

% Set up parameters for the *early 2000s* case

% Parameter Values of oil supply shocks, early 2000s
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p =  0.0655*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8158;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks

% Parameter Values of TFP shocks *early 2000s*
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 0.2707*0.3694*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8466;   % serial correlation of transitory TFP shocks
sigz_t = 0.3694*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep  = []; % not useful here

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *early 2000s parameters - with inventories* solution:')
disp('Will store std deviation(GDP) and std (spot prices) to construct the table :')
simule


table5panelB(1,1) = mvar(iy);
table5panelB(2,1) = mvar(ispot);
table5panelB

% Now working on the *no inventories* case for the same shock process
% only need to change "answinv"
answinv = 0;
% answinv = 1 implies solution where oil inventoires are present

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *early 2000s parameters - no inventories* solution:')
disp('Will select std deviation(GDP) and std (spot prices) and express the previous ones as a ratio of the these to construct the table :')
simule


table5panelB(1,1) = table5panelB(1,1)/mvar(iy);
table5panelB(2,1) = table5panelB(2,1)/mvar(ispot);
table5panelB


% Now on to the second column of Table 5 Panel B : *late 2000s* shocks 


% Set up parameters for the *late 2000ss* case

% Parameter Values of oil supply shocks, *late 2000s* values
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p = 0.2408*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.8294;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;     % std. deviation of innovations to those shocks
% Parameter Values of TFP shocks, *late 2000s* values
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 2.4926*0.2592*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.8463;   % serial correlation of transitory TFP shocks
sigz_t = 0.2592*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep  = []; % not useful here

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *late 2000s parameters - with inventories* solution:')
disp('Will store std deviation(GDP) and std (spot prices) to construct the table :')
simule


table5panelB(1,2) = mvar(iy);
table5panelB(2,2) = mvar(ispot);
table5panelB

% Now working on the *no inventories* case for the same shock process
% only need to change "answinv"
answinv = 0;
% answinv = 1 implies solution where oil inventoires are present

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *early 2000s parameters - no inventories* solution:')
disp('Will select std deviation(GDP) and std (spot prices) and express the previous ones as a ratio of the these to construct the table :')
simule


table5panelB(1,2) = table5panelB(1,2)/mvar(iy);
table5panelB(2,2) = table5panelB(2,2)/mvar(ispot);
table5panelB


% Now on to the third column of Table 5- Panel B : *mid 2010s* shocks 


% Set up parameters for the *mid 2010s* case

% Parameter Values of oil supply shocks, "mid 2010s" values
rhoe_p = 0.99999;   % serial correlation of persistent oil-supply shocks
sige_p = 0.1344*0.01; % std. deviation of innovations to those shocks
rhoe_t = 0.6899;   % serial correlation of transitory oil-supply shocks
sige_t = 0.01;    % std. deviation of innovations to those shocks
% Parameter Values of TFP shocks, "mid 2010s" values
rhoz_p = 0.99999;   % serial correlation of persistent TFP shocks
sigz_p = 0.5624*0.6159*0.01; % std. deviation of innovations to those shocks
rhoz_t = 0.9322;   % serial correlation of transitory TFP shocks
sigz_t = 0.6159*0.01;    % std. deviation of innovations to those shocks

whatinfo = 2;
% whatinfo = 1 implies "complete information" solution as defined in the paper; whatinfo = 2 implies "learning" solution 
answinv = 1;
% answinv = 1 implies solution where oil inventoires are present
answlearn = 1;
% answlearn = 1 indicates the usual learning solution. "answlearn" = 0 implies the learning solution used only to construct Table 5 

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];

gainkeep  = []; % not useful here

% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *mid 2010s parameters - with inventories* solution:')
disp('Will store std deviation(GDP) and std (spot prices) to construct the table :')
simule


table5panelB(1,3) = mvar(iy);
table5panelB(2,3) = mvar(ispot);
table5panelB

% Now working on the *no inventories* case for the same shock process
% only need to change "answinv"
answinv = 0;
% answinv = 1 implies solution where oil inventoires are present

param = [rhoe_p sige_p rhoe_t sige_t rhoz_p sigz_p rhoz_t sigz_t answinv whatinfo answlearn];


% Call the program that computes the s.s. and the linearization

lin
system = 'lin';

% Call the KW(2002) program

[M,PI,G]=resolkw(system,track,comp,param,gainkeep);

disp('Working on Panel B of the table:')
disp('Simulating 500 data samples for the *mid 2010s parameters - no inventories* solution:')
disp('Will select std deviation(GDP) and std (spot prices) and express the previous ones as a ratio of the these to construct the table :')
simule


table5panelB(1,3) = table5panelB(1,3)/mvar(iy);
table5panelB(2,3) = table5panelB(2,3)/mvar(ispot);
table5panelB

disp(' ********  ')
disp(' All Done  ')
disp('           ')
disp(' Here is Panel A of the table : ')
disp(table5panelA)
disp('           ')
disp(' Here is Panel B of the table : ')
disp(table5panelB)




